In [1]:
%pylab inline
import matplotlib.pyplot as plt

#import collections
Populating the interactive namespace from numpy and matplotlib

Phonon spectrum for cubic systems including first- and second-neighbor interactions in the potential between atoms

Jenő Sólyom: Fundamentals of the Physics of Solids, Volume 1: Structure and Dynamics (A modern szilárdtestfizika alpajai I. A szilárd testek szerkezete és dinamikája)

See: page 350. Eq. (11.2.56) for simple cubic structure including 1st and 2nd neighbor interactions

In [12]:
# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 12  #    12 nagy abrahoz
yfig_meret= 7    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 20  #  30 nagy abrahoz
legend_meret= 21   #  30 nagy abrahoz
In [3]:
def Dm_op(qv,K1,K2,delta1,delta2):
    
    ## q-points are in units of 2 pi/a, qv

    Dm1= array(zeros((3,3)))
    for d in delta1:
        Dm1 = Dm1 + outer(d,d)/dot(d,d)*(1-exp(1j*2*pi*dot(qv,d)))
    
    Dm2= array(zeros((3,3)))
    for d in delta2:
        Dm2 = Dm2 + outer(d,d)/dot(d,d)*(1-exp(1j*2*pi*dot(qv,d)))
    
    Dm= K1*Dm1 + K2*Dm2
    return  Dm
In [4]:
#  simple cubic 

qv = [1.2,1.2,0]
K1=2
K2=1

deltaSC1= array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]), 
                 array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])

deltaSC2=array([array([0, 1, 1]), array([0, 1, -1]), array([0, -1, 1]), 
                array([0, -1, -1]), array([1, 0, 1]), array([1, 0, -1]), 
                array([1, 1, 0]), array([1, -1, 0]), array([-1, 0, 1]), 
                array([-1, 0, -1]), array([-1, 1, 0]), array([-1, -1, 0])])

real(Dm_op(qv,K1,K2,deltaSC1,deltaSC2));

Phonon spectrum for simple cubic lattice

Brillouin zone with high symmetry $\mathbf{k}$ points

Dispersion relations for the lattice vibrations along the four characteristic directions of the Brillouin zone in a simple cubic crystal with a monatomic basis

see page 354, Fig. 11.10.

In [13]:
K1=1
K2=0.5

Npoints = 50;  ## hany pontbol keszul a gorbe
ymax = 3.0;

## k-points are in units of 2 pi/a

## high symmetry points in the Brillouin zone
## FONTOS: a fenti abra szerinti pontokat kell venni,
## a szomszedos spec. k pontokat kell venni.
Gp = array([0,0,0])
Xp = 0.5*array([0,1,0])
Rp = 0.5*array([1,1,1])
Mp = 0.5*array([1,1,0])

##  G-X-M-G-R line
specpoints = array([Gp,Xp,Mp,Gp,Rp])

specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Mp):"$M$",
               tuple(Rp):"$R$"}

klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
    sl = specpoints[i+1]-specpoints[i]  
    sl_hossz = sqrt(dot(sl,sl))
    dk = sl/Npoints
    for num in arange(0,Npoints):
        # k vector along the line given by speclines:
        klist.append(specpoints[i] + num*dk)
        # length of the k vector and shifted: 
        tmp +=  sqrt(dot(dk,dk)) 
        kk.append(tmp) 
    specpos.append(tmp)

klist.append(klist[-1]+dk)   

enlist = []
for kv in klist:
    enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaSC1,deltaSC2))))
    
enlist= array(enlist)

eigszam= len(enlist[0,:])
  
specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]

    
# drawing the figure
figsize(xfig_meret,yfig_meret)

for ii in range(0,eigszam):
    plot(kk,enlist[:,ii],color='r')

i = 0
for xc in specpos:
    axvline(x=xc,color='b')
    text(xc,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
    i=i+1   
    
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
xticks([], [])

xlim(0,specpos[-1])
ylim(0,ymax)

cim = "phonon for sc, $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)

grid();
    
#savefig('Fig_phonon_simple_cubic_1-2_neighbors.pdf');  # Abra kimentese

Brillouin zone with high symmetry $\mathbf{q}$ points for a fcc lattice

see

Dispersion relation for a face-centered cubic gold crystal

Jenő Sólyom: Fundamentals of the Physics of Solids, Volume 1: Structure and Dynamics (A modern szilárdtestfizika alpajai I. A szilárd testek szerkezete és dinamikája) see here on page 355. Fig. 11.11.

In [17]:
## fcc lattice:

K1=5
K2=0.5

Npoints = 500;  ## hany pontbol keszul a gorbe
ymax = 3.5;

## k-points are in units of 2 pi/a

## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni 
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok, 
## es igy jo az abra)
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])

##  G-X-W-L-G-K-U-X line
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Up,Xp])
#specpoints =  array([Lp,Gp,Xp,Up,Kp,Gp,Wp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp])  #  Solyom konyvben: II. kotet, 20.1. abra
                                         #  Az  U,K pontok egybe!!!!! 

specpoints = array([Gp,Xp,Wp,Xp,Up,Kp,Gp,Lp])

## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy 
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!

###  kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
               tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}

klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
    if not i==4:   ########## Az  U,K pontok egybe!!!!! 
        # a i== 4 jeloli, hogy hol van az Up pont, 
        # azaz az Up indexe a 'specpoints' listaban. 
          
        sl = specpoints[i+1]-specpoints[i]  
        sl_hossz = sqrt(dot(sl,sl))
        dk = sl/Npoints
        for num in arange(0,Npoints):
           # k vector along the line given by speclines:
             klist.append(specpoints[i] + num*dk)
           # length of the k vector and shifted: 
             tmp +=  sqrt(dot(dk,dk)) 
             kk.append(tmp) 
        specpos.append(tmp)

klist.append(klist[-1]+dk)   

deltaFCC1= array([array([0, 1/2, 1/2]), array([0, 1/2, -1/2]), array([0, -1/2, 1/2]), 
                  array([0, -1/2, -1/2]),
                  array([1/2, 0, 1/2]), array([1/2, 0, -1/2]), array([-1/2, 0, 1/2]), 
                  array([-1/2, 0, -1/2]),
                  array([1/2, 1/2, 0]), array([1/2, -1/2, 0]), array([-1/2, 1/2, 0]), 
                  array([-1/2, -1/2, 0])])

deltaFCC2=array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]), 
                 array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])

enlist = []
for kv in klist:
    enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaFCC1,deltaFCC2))))
    
enlist= array(enlist)

eigszam= len(enlist[0,:])
  
specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
  
# drawing the figure
figsize(xfig_meret,yfig_meret)

for ii in range(0,eigszam):
    plot(kk,enlist[:,ii],color='r')

specNev = ["$\Gamma$","$X$","$W$","$X$","$K$","$\Gamma$","$L$" ]
#specpoints = array([Gp,Xp,Wp,Xp,Up,Kp,Gp,Lp])

               
i = 0
for xc in specpos:
    axvline(x=xc,color='b')
    #text(xc,-0.5,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
    text(xc,-0.5,specNev[i],fontsize=xylabel_meret)
    i=i+1   
    
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
xticks([], [])

xlim(0,specpos[-1])
#ylim(0,ymax)

cim = "phonon for fcc (Au) $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)

grid();


#savefig('Fig_phonon_fcc_1-2_neighbors.pdf');  # Abra kimentese

Phonon spectrum for bcc lattice

Brillouin zone with high symmetry $\mathbf{k}$ points

In [16]:
# bcc lattice 

K1=1
K2=0.5

deltaBCC1 = array([array([1/2, 1/2, 1/2]), array([1/2, 1/2, -1/2]), 
                   array([1/2, -1/2, 1/2]), array([1/2, -1/2, -1/2]),
                   array([-1/2, 1/2, 1/2]), array([-1/2, 1/2, -1/2]), 
                   array([-1/2, -1/2, 1/2]), array([-1/2, -1/2, -1/2])])

deltaBCC2 = array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]), 
                 array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])

Npoints = 50;  ## hany pontbol keszul a gorbe
ymax = 3.0;

## k-points are in units of 2 pi/a

## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni 
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok, 
## es igy jo az abra)
Gp = array([0,0,0])
Hp = 0.5*array([0,1,0])
Pp = 0.25*array([1,1,1])
Np = 0.25*array([1,1,0])

##  G-H-N-G-P-H line
specpoints = array([Gp,Hp,Np,Gp,Pp,Hp])


specNev_map = {tuple(Gp):"$\Gamma$",tuple(Hp):"$H$",tuple(Pp):"$P$",
               tuple(Np):"$N$"}
            
figsize(xfig_meret,yfig_meret)

klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
    sl = specpoints[i+1]-specpoints[i]  
    sl_hossz = sqrt(dot(sl,sl))
    dk = sl/Npoints
    for num in arange(0,Npoints):
        # k vector along the line given by speclines:
        klist.append(specpoints[i] + num*dk)
        # length of the k vector and shifted: 
        tmp +=  sqrt(dot(dk,dk)) 
        kk.append(tmp) 
    specpos.append(tmp)

klist.append(klist[-1]+dk)   

enlist = []
for kv in klist:
    enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaBCC1,deltaBCC2))))
    
enlist= array(enlist)

eigszam= len(enlist[0,:])
  
#specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]

    
# drawing the figure
figsize(xfig_meret,yfig_meret)

for ii in range(0,eigszam):
    plot(kk,enlist[:,ii],color='r')

i = 0
for xc in specpos:
    axvline(x=xc,color='b')
    text(xc,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
    i=i+1   
    
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
xticks([], [])

xlim(0,specpos[-1])
#ylim(0,ymax)

cim = "phonon for bcc, $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)

grid();
    
#savefig('Fig_phonon_bcc_1-2_neighbors.pdf');  # Abra kimentese
In [ ]: